The earth's electromagnetic response is commonly measured both offshore and on land in an effort to identify possible areas of hydrocarbon accumulation. Controlled-source and/or magnetotelluric data are used to derive a resistivity model of the subsurface, (Ellingsrud et al., 2002; Srnka et al., 2006) which is used to identify the potential areas of hydrocarbon accumulation in the subsurface. Such subsurface models can be equivalently described in terms of their conductivities.
Measured electromagnetic data have a complicated, non-linear dependence on the subsurface resistivity. Subsurface resistivities are typically derived from electromagnetic data by iterating a linearized inversion designed to minimize an objective function that characterizes the error or misfit between acquired and synthesized data (Newman and Alumbaugh, 1997 and Carazzone et al., 2005). So, in FIG. 1, a typical objective function might be
                    ϕ        =                                            1              2                        ⁢                                          ∑                i                            ⁢                                                                    (                                                                  W                        i                                            ⁡                                              (                                                                              d                            i                            meas                                                    -                                                                                    d                              i                              synth                                                        ⁡                                                          (                                                              m                                →                                                            )                                                                                                      )                                                              )                                    T                                ⁢                                  (                                                            W                      i                                        ⁡                                          (                                                                        d                          i                          meas                                                -                                                                              d                            i                            synth                                                    ⁡                                                      (                                                          m                              →                                                        )                                                                                              )                                                        )                                                              +                                    1              2                        ⁢            λ            ⁢                                                  ⁢                                          F                T                            ⁡                              (                                  m                  →                                )                                      ⁢                          F              ⁡                              (                                  m                  →                                )                                                                        (        1        )            Here, the di are measured or synthetic electromagnetic field values and the sum is over all data, i, including locations, vector components, and frequencies. The Wi are data weights chosen to compensate the data for varying magnitudes and uncertainties and {right arrow over (m)} is the vector of unknown subsurface parameters such as conductivity, log conductivity, or resistivity. The vector, {right arrow over (m)}, may further include certain unknown or poorly-determined parameters of the survey, such as the orientation of receiver antennae on the seafloor. The superscript T represents the matrix transpose. The function, F, is chosen to steer the inversion toward an a priori notion of the subsurface parameters. The trade-off parameter, λ, controls the relative strength of the data mismatch and a priori terms and is typically reduced as the data mismatch decreases during the inversion. As shall be seen below, it is possible to introduce additional terms into equation (1), with additional trade-off parameters, λi, to encourage the inversion toward models that match several a priori criteria.
The data mismatch term in the objective function of equation (1) is usually made quadratic in {right arrow over (m)} by assuming that the model change from one update in an iterative inversion process is small. To make the total objective function, φ, quadratic in {right arrow over (m)} and thereby simplify the least-squares minimization, the regularization function F is typically chosen to be linear in {right arrow over (m)}. A linear regularization function can be expressed in matrix notation as,F({right arrow over (m)})=C·{right arrow over (m)}  (2)The matrix C represents the linear operation of the regularization on the inversion parameters, {right arrow over (m)}.
The trade-off parameter, λ, controls the relative strength of the a priori term and the data mismatch and is typically reduced as the data mismatch decreases during the inversion. λ will most typically be adjusted during the inversion iterations so that the regularization term in the objective function is about one order of magnitude smaller than the data mismatch term and preferably between ⅕th and 1/100th of the data mismatch term.
The efficient implementation of a nonlinear conjugate gradient method on massively-parallel computers and, in particular, the separation of the grids used to compute the electromagnetic fields from the grids used to represent the subsurface model have been discussed by Commer and Newman (2008) and by Newman and Alumbaugh (1997, 2000). In addition to the conjugate gradient technique, various other inversion methods have been used by geophysicists to invert either electromagnetic or seismic data. Such methods include differential semblance optimization (Symes and Carazzone, 1991), genetic algorithms, (Stoffa and Sen, 1991), simulated annealing, (Landa, et al, 1989), maximum smoothness (Constable, et al, 1987), and Gauss-Newton (Hu et al., 2007).
Non-uniqueness is a feature common to most geophysical data inversions. This is especially true for controlled source electromagnetic (“CSEM”) inversion given the exponential decay of the data, the low frequencies involved, and the limited data coverage. The regularization term in equation (1) is used to mitigate the non-uniqueness in the solution. For the case of linear regularization function, the operator C in equation (2) is normally set to be either first-order or second-order spatial derivatives (Newman and Alumbaugh (1997, 2000); Constable, et al, 1987), which leads to a spatially smooth solution. If, instead, the term is implemented as,F({right arrow over (m)})C0·({right arrow over (m)}−{right arrow over (m)})  (3),where C0 is a diagonal matrix and {right arrow over (m)}0 is an a priori estimate of the model parameters, then the inversion will be prejudiced toward models that are in the neighborhood of the a priori estimate {right arrow over (m)}0.
Non-linear regularizations have also been used in geophysical inversion. One such example is (Portniaguine, 1999; Blaschek, 2008),
                                                        F              T                        ⁡                          (                              m                →                            )                                ⁢                      F            ⁡                          (                              m                →                            )                                      =                                                                              (                                      ∇                                                                                  ⁢                                          m                      →                                                        )                                T                            ⁢                              (                                  Δ                  ⁢                                                                          ⁢                                      m                    →                                                  )                                                                                                          (                                          ∇                                                                                          ⁢                                              m                        →                                                              )                                    T                                ⁢                                  (                                      Δ                    ⁢                                                                                  ⁢                                          m                      →                                                        )                                            +                              β                2                                              .                                    (        4        )            This regularization is designed to seek solutions with a minimum number of sharp boundaries by adjusting the value of parameter β.
Even with these different types of regularization, the resistivity images from CSEM data inversions generally still lack structural consistency with subsurface structure derived from other sources, such as seismic data. While the electrical and elastic properties of the earth may differ in detail at the reservoir scale, they will generally be correlated over larger distances since they arise from the same subsurface region, depositional processes, source rocks, and pore fluids. This structure inconsistency makes geological interpretation difficult and the prediction of hydrocarbon presence less reliable.
Another important aspect of subsurface imaging is rock physical property anisotropy. The subsurface resistivity may be anisotropic. That is, the electrical resistivity of the earth may depend on the direction of current flow (Lu and Xia, 2007; Jupp and Vozoff, 1997). Likewise, seismic data are also acquired and imaged to form a model of the earth's elastic properties, which may also be anisotropic. Since the subsurface rock properties are both anisotropic and inhomogeneous over the distances probed by geophysical data, it is often difficult to unravel intrinsic from structurally-induced anisotropy. The intrinsic rock-property anisotropy represents the anisotropy of the rock itself, which exists normally at micro scale much smaller than that probed by seismic or CSEM data. A typical example is anisotropy due to the microscopic, plate-like structures in shales and other clay-rich minerals. This intrinsic rock anisotropy may and may not show up at large scale depending how the rock fine grains are aligned in space. So, shales deposited gently in large sedimentary beds over long periods of time may manifest their intrinsic anisotropy in electromagnetic data while sand deposited in a more rapid and chaotic fashion may appear isotropic in electromagnetic data. Structurally-induced anisotropy represents the effective anisotropy at the scale probed by a particular geophysical data type due to the rock property variation in space within the probed scale. Normally, the subsurface shows structurally-induced anisotropy at the probing scale due to the inhomogeneous rock property at the subsurface even the rocks are intrinsically isotropic at fine scale. Thus, alternating beds of sandstone and limestone, while intrinsically isotropic, may manifest themselves as an anisotropic composite in electromagnetic data. While the symmetry axes of the electrical properties may be oriented differently than those of the elastic properties, these orientations are intrinsically correlated. In principle, the full anisotropic resistivity, including both the orientations of the symmetry axes and the principle resistivities (resistivities along the symmetry axes) can be treated as unknown parameters in an anisotropic inversion (Pain, et al, 2003). Such straightforward electrical inversion can fail to predict the orientations of the symmetry axes (even on synthetic data from a simple model) when the orientations vary in space (Pain, et al, 2003). It is difficult to derive reliable orientations of symmetry axes and principle resistivities from the measured CSEM data in geological settings encountered in hydrocarbon exploration, and it is even more difficult to correlate the derived orientations to that from other sources, such as seismic data.
In any event, the problems of inverting CSEM data to a resistivity model that is consistent with the subsurface structure derived from other sources, such as seismic data, can not be achieved by those common regularizations without utilizing the a priori structure information.
In order to understand the differences in detail and to render a useful geological interpretation from the resistivity and elastic models, it is very valuable to have methods that will invert CSEM data in a way that is globally consistent with seismic or other geophysical measurements, and distinguish the intrinsic subsurface anisotropy from the effects of geological structures and inhomogeneity.
One direct solution that has been proposed for obtaining a resistivity image consistent with another property, such as the seismic velocity, is to link the electrical rock property directly to said other property (such as elastic properties) by using some relationships (Archie, 1942; Xu and White, 1995), and to perform joint inversion of the CSEM data and geophysical data that measure said other property. For the case of joint electromagnetic and seismic inversion (Hoversten, et al 2006), the inversion parameter vector, {right arrow over (m)}, includes parameters from which both resistivities and subsurface elastic properties such as the density and seismic velocities can be derived by using known relationships. The data vectors for joint inversion, di, include both electromagnetic field measurements and seismic measurements such as pressure, velocity, or acceleration. The objective function measures both the misfit between synthetic and measured electromagnetic data as well as the misfit between synthetic and measured seismic data; and, the function F in equation (1) is chosen to enforce the smoothness requirements by means of a simple, first-order gradient of the parameters. As the electrical and elastic properties of the subsurface are updated to reduce the sum of the data misfit and model smoothness terms, the updated parameters are additionally constrained to satisfy rock-physics relationships. One example of these additional constraints would be to require that the electrical and elastic properties be derivable from unique values for porosity, water saturation, and shale volume. The seismic and resistivity images are thereby linked by the rock physics relationships during this type of inversion.
Due to the complexity of the subsurface rock, it is very challenging to build rock physics relationships that work for the wide variety of rock mixtures. Chen and Dickens (Chen and Dickens, 2007) have shown that errors in the rock physics relationships have a big impact on the reliability of the inversion results. Moreover, different geophysical data types respond to changes in subsurface properties over very different length scales, making it ambiguous how the properties should be averaged into the rock-property relationships. These errors and complexities make the practical application of this direct method of joint inversion difficult.
The standard method (see, for example, A. Lovatini, K. Umbach, and S. Patmore, “3D CSEM inversion in a frontier basin offshore west Greenland,” First Break 27, 95-98 (2009)) for inverting electromagnetic data to resistivity images that are consistent with models of other rock properties (such as a seismic reflectivity image, seismic velocity, gravity model, magnetic model, or well logs) is to begin the electromagnetic data inversion with a resistivity model that is already at least structurally consistent with the non-electromagnetic data. The inversion usually breaks this consistency due to the discrepancy between the starting model and the actual subsurface resistivity, limited data coverage, resolution, and non-uniqueness. The inversion result is then re-interpreted and modified in light of the seismic images. The modified resistivity model then serves as the starting model for a subsequent inversion. If successful, the method converges to a resistivity model that is consistent with both the seismic image and the electromagnetic data. This method is both labor-intensive and computationally inefficient. It becomes impractical or can fail to converge to the correct answer when more than two types of data are included.
In an alternate approach to the rock-physics-relationship based joint inversion, Gallardo et al. (2004, 2005) introduced the joint inversion of zero-frequency (i.e., direct current) electrical measurements and seismic refraction travel time measurements with cross-gradient constraints. In addition to the usual regularization (smoothing) terms, the side condition or cross-gradient constraint,∇{right arrow over (m)}e×∇{right arrow over (m)}s=0  (5)is applied to the electrical parameter vector {right arrow over (m)}e and the seismic parameter vector {right arrow over (m)}s. This cross-product condition is designed to ensure the geometrical similarity between the inverted results of {right arrow over (m)}e and {right arrow over (m)}s.
As a refinement to the method of Gallardo et al., Hu et al. (2007) separate the joint inversion of electric and seismic data into cyclic seismic and electromagnetic inversion steps. After the electromagnetic parameters have been updated, the seismic parameters are inverted and updated, followed by an inversion and update of the electromagnetic parameters and so on. Each inversion is based on an objective function which includes a cross-gradient-type a priori term from the parameters corresponding to other data type
                              ϕ          e                =                                            1              2                        ⁢                                          ∑                i                            ⁢                                                (                                                            W                      i                                        ⁡                                          (                                                                        d                          i                          meas                                                -                                                                              d                            i                            synth                                                    ⁡                                                      (                                                                                          m                                →                                                            e                                                        )                                                                                              )                                                        )                                2                                              +                                    1              2                        ⁢                                                            λ                  1                                ⁡                                  [                                      C                    ·                                          (                                                                                                    m                            →                                                    e                                                -                                                                              m                            →                                                                                e                            ,                            0                                                                                              )                                                        ]                                            2                                +                                    λ              2                        ⁢                                                                                                ∇                                                                                  ⁢                                                                  m                        →                                            e                                                        ×                                      ∇                                                                                  ⁢                                                                  m                        →                                            s                                                                                                  2                                                          (                  6          ⁢                                          ⁢          a                )                                          ϕ          s                =                                            1              2                        ⁢                                          ∑                i                            ⁢                                                (                                                            W                      i                                        ⁡                                          (                                                                        d                          i                          meas                                                -                                                                              d                            i                            synth                                                    ⁡                                                      (                                                                                          m                                →                                                            s                                                        )                                                                                              )                                                        )                                2                                              +                                    1              2                        ⁢                                                            λ                  1                                ⁡                                  [                                      C                    ·                                          (                                                                                                    m                            →                                                    s                                                -                                                                              m                            →                                                                                s                            ,                            0                                                                                              )                                                        ]                                            2                                +                                    λ              2                        ⁢                                                                                                                        ∇                                                                                          ⁢                                                                        m                          →                                                s                                                              ×                                          ∇                                                                                          ⁢                                                                        m                          →                                                e                                                                                                              2                            .                                                          (                  6          ⁢                                          ⁢          b                )            This approach allows Hu et al. to apply the more sophisticated Gauss-Newton method (Hu et al., 2007) separately for each data type, updating the electrical (or seismic) parameters in the seismic (or electrical) objective function after each iteration of the electrical (or seismic) inversion.
In practical application, it is very difficult for the two images to converge to the same structure due to factors such as the limited data coverage, rapid decay of the electric fields, subsurface anisotropy, unknown orientations of symmetry axes, and errors in the initial electrical and seismic models. These factors can lead to model gradients that are very poorly behaved (poorly correlated with the seismic data) during the iterations. Also, the subsurface anisotropy, which is important for the inversion of electromagnetic data, is not addressed by this cross-gradient method.
Along the line of adding structure information in the objective function, a method of seismically-constrained resistivity inversion has been introduced by Saunders et al, (Saunders et al, 2005) for the inversion of zero-frequency (i.e., direct current (DC)) electrical measurements. This work follows earlier publications on zero-frequency (DC) resistivity inversion (Kaipio et al, 1999; Herwanger et al, 2002; Pain et al, 2003). The a priori structure is derived from a known seismic velocity field and is incorporated into a constraint term in the objective function in the form of½(∇{right arrow over (m)}T)(∇{right arrow over (m)})  (7)The idea is to apply small or weak regularization in the direction and area where there is large curvature in the structure field, thus permitting larger changes to the electrical properties. The symmetric tensor  in equation (7) is obtained from the curvature of the known structure field, such as a seismic velocity field. More specifically, , is constructed from the eigenvalues and eigenvectors of the normalized Hessian of the known structure field, which, in turn, involves all the second-order spatial derivatives of the structure field.
The methods of Saunders et al, Paipio et al, Herwanger et al, and Pain et al focus on zero-frequency (i.e., direct current) electrical measurements. The full anisotropic resistivity, (including both the principle resistivity values and the orientations of the symmetry axes) were assumed to be derivable from the inversion of DC electrical measurements, and orientations of the symmetry axes were assumed to follow the same known structure. These authors, moreover, specifically add a regularization term to penalize the intrinsic anisotropy of the unknown material, apparently on the theory that allowing the intrinsic resistivities to vary freely can result in erroneous solutions and that the most realistic solutions are those with the smallest anisotropy. Based on the present inventors' experience with both surface electromagnetic and well-log data, this approach is rejected, and a term to penalize the intrinsic anisotropy is not added in the present invention.
Since the subsurface rock properties are both anisotropic and inhomogeneous over the distances probed by geophysical data, it is often difficult to unravel intrinsic anisotropy (a property of the rocks) from structurally-induced anisotropy (an aggregate property of rocks). Even though the orientations of the symmetry axes of electrical properties may be intrinsically correlated to those of the elastic properties, they may differ when averaged over the length scales appropriate to different geophysical data types. For example, the intrinsic anisotropy of shale beds may have been tilted relative to the horizontal during lithification, so that they manifest tilted anisotropy when measured by well-logging tools in a vertical wellbore. Meanwhile, adjacent sand beds would likely be isotropic. Surface electromagnetic data, operating at longer wavelengths than well logs, will sense a mixture of the tilted intrinsic anisotropy from the shales and a structural anisotropy from the interbedding of shales and sands. Simply forcing both the resistivities and orientations of the electrical symmetry axes to follow the same a priori structure during CSEM inversion may not generate geologically-consistent and reliable orientations of the electrical symmetry axes. The error in symmetry axis orientations will translate to error in the resistivity images, which can impact the interpretation of the subsurface. The non-uniqueness problem exists in most geophysical data inversions, and is especially true for CSEM data due to the low resolution and exponential decay of the fields. All these factors combined with the complex geological settings encountered in hydrocarbon exploration make the prediction of the full anisotropic resistivities, especially the orientations of the symmetry axes, difficult using CSEM data alone. It is therefore desirable to have an effective and practical method for inverting CSEM data to a resistivity model that is consistent with a subsurface structure derived from other sources, such as seismic data.